################################################################################
#                                                                              #
#        Replication script for analysis contained in:                         #
#        If the Jury Only Knew: The Effect of Omitted Mitigation               #
#        Evidence on the Probability of a Death Sentence,                      #
#        By Barry C. Edwards, J.D., Ph.D.                                      #
#                                                                              #
#        Published in Virginia Journal of Social Policy & the Law              #
#        Spring 2025                                                           #
#                                                                              #
################################################################################


################## Set-up machine and load datasets ############################
rm(list=ls())  # starting with clean environment

# This script requires the sate and RCPA3 packages.
# The next lines installed them if you don't have them.
if(require("sate", quietly = TRUE)==FALSE) install.packages("sate")
if(require("RCPA3", quietly = TRUE)==FALSE) install.packages("RCPA3")
# make sure sure you're using sate version 2.3.0 or higher
if(packageVersion("sate") < "2.3.0") install.packages("sate")

# load the installed sate and RCPA3 packages.
library(sate)  
library(RCPA3)

# Load dataset. It is saved as "observed.deliberations" in sate package.
all_juries_data <- observed.deliberations
 

########################## Set-up for analysis #################################

all_juries_data$lenient_verdict <- 1 - all_juries_data$guilty_verdict

# the comparison_table is not reported, the results are is displayed in figure 1
all_juries_data$initial_lenient_votes <- all_juries_data$jury_size - round(all_juries_data$prop_jurors_g * all_juries_data$jury_size)
all_juries_data$initial_lenient_fraction <- all_juries_data$initial_lenient_votes / all_juries_data$jury_size

comparison_table <- RCPA3::compmeansC(dv=all_juries_data$lenient_verdict, 
                                      iv=factor(all_juries_data$initial_lenient_fraction), 
                                      plot=F)
observed_life_probs <- as.numeric(comparison_table[1:13,"Mean"])
# rows 1:13 correspond to 0: 12 votes for death sentence


########################## Figure 1 #################################

# Figure 1 is bar plot of proportion of lenient verdicts by fraction for leniency
# I have commented out the png() and dev.off() commands
# Uncomment these commands to generate .png version of Figure 1.

png(file="basic observed bar plot for VJSPL.png", width=4, height=3, 
    units="in", type="cairo", pointsize=11, res=300, antialias="default")

par(mar=c(3, 3.25, 0, 0), family="serif", mfrow=c(1,1), omi=base::c(0.1,0.1,0.1,0.1), tck=-.02)
yaxis_labels <- c("0.0", "0.1", "0.2", "0.3", "0.4", "0.5", "0.6", "0.7", "0.8", "0.9", "1.0")
fraction_labels <- c(0, expression(frac(1,12)), expression(frac(1, 6)),  
                     expression(frac(1, 4)), expression(frac(1, 3)), expression(frac(5, 12)),  
                     expression(frac(1, 2)), expression(frac(7, 12)), expression(frac(2, 3)),  
                     expression(frac(3, 4)), expression(frac(5, 6)), expression(frac(11, 12)), 1)
barplot(height = observed_life_probs, cex.names=.8, ann = T, 
        col = c("gray90"), 
        args.legend=list(x=7.5, y=1, cex=.9, box.col = "white", box.lty = 0),
        xlab="", ylab="", axes=F, ylim=c(0,1.05), space = .25)
axis(side=2, at=seq(0, 1, by=.10), labels = yaxis_labels, las=2, cex.axis=.8)
mtext(side=1, at=.75 + 0:12*1.25, text = fraction_labels, cex=c(.8,rep(.6,11),.8), line = .25, padj = .5)
mtext(side = 1, text = "Initial Fraction of Jurors for Leniency", cex = .9, line = 2)
mtext(side = 2, text = "Proportion of Lenient Verdicts", cex = .9, line = 2)

# dev.off()


############## Table 3: Results of logistic regression analysis ##################

# For logistic regression analysis, there are more than 13 possible fractions
# due to averaging some reported first poll results from a field study
all_juries_data$initial_lenient_fraction <- 1 - all_juries_data$prop_jurors_g

model1 <- logregC(lenient_verdict ~ initial_lenient_fraction, data=all_juries_data, 
                  fit.stats = T, pre = T)
model2 <- logregC(lenient_verdict ~ initial_lenient_fraction + death_penalty + six_person_jury, 
                  data=all_juries_data, fit.stats = T, pre = T) 


######## Figure 2: Line chart of predicted probability by verdict fraction ######

# based on predicted probabilities from logistic regression analysis, Model 1
predicted_probs = predict.glm(model1, newdata=data.frame(initial_lenient_fraction=0:12/12), 
                              type="response", se.fit=T)

# I have commented out the png() and dev.off() commands
# Uncomment these commands to generate .png version of Figure 2.
# png(file="predicted probabilities plot for VJSPL.png", width=4, height=3, 
#     units="in", type="cairo", pointsize=11, res=300, antialias="default")

yaxis_labels <- c("0.0", "0.1", "0.2", "0.3", "0.4", "0.5", "0.6", "0.7", "0.8", "0.9", "1.0")
par(mar=c(3, 3.25, 0, 0), family="serif", mfrow=c(1,1), omi=base::c(0.1,0.1,0.1,0.1), tck=-.02)
plot(x="", y="", xlim=c(0,12), xlab="",
     ylab="", axes=F, ylim=c(0,1.05))
axis(side=2, at=seq(0, 1, by=.10), labels = c(yaxis_labels), las=2, cex.axis=.8)
axis(side=1, at=seq(0, 12, by=1), labels = rep("", 13), cex.axis=.8)
mtext(side=1, at=0:12, text = fraction_labels, cex=c(.8,rep(.6,11),.8), line = 0.4, padj = .5)
mtext(side = 1, text = "Initial Fraction of Jurors for Leniency", cex = .9, line = 2)
mtext(side = 2, text = "Probability of a Lenient Verdict", cex = .9, line = 2)
box(bty="L")
lines(x=0:12, y=predicted_probs$fit)
points(x=0:12, y=predicted_probs$fit, pch=21, col="black", bg="gray")

# dev.off()


